Change Detection - Auswertung der Klassifikationsdaten

Die Auswertung

#
# Type: script
# Name: change_detection_workflow.R
# Description:  basic reproducible Change Detection workflow for Sentinel data
#               data download -> digitizing training areas ->
#               extracting and preparing training data ->
#               example classifications: cluster analysis, random forest, MaxLike
#               quality assesment and comparison with basic tools, kappa, and motif
#               basic visualisation examples with tmap and mapview
# Dependencies:
# Output: original sentinel tile
#         AOI window of this tile (research_area)
#         training areas
# Copyright: GPL (>= 3), Chris Reudenbach, creuden@gmail.com
#

#--- laden der notwendigen Bibliotheken
# Achtung Pakete müssen evtl. manuell installiert werden
library(envimaR)
library(tmaptools)
library(motif)
library(kableExtra)
library(tmap)
library(stars)

#--- Schalter für den Download der sentinel daten
get_sen = FALSE

#--- schalter ob digitalisiert werden muss falls er auf FALSE gesetzt ist werden die
# (zuvor erstellten und gesicherten Daten ) im else Teil der Verzweigung eingelesen
digitize = FALSE

## setzen des aktuellen Projektverzeichnisses (erstellt mit envimaR) als root_folder
#root_folder = find_rstudio_root_file()
root_folder = "~/edu/geoinfo/"

# einlesen des zuvor erstellten Setup-Skripts
source(file.path(root_folder, "src/functions/000_setup.R"))
source(file.path(root_folder, "src/functions/fun_helper.R"))
source(file.path(root_folder, "src/functions/fun_kappa.R"))
source(file.path(root_folder, "src/functions/fun_tmap.R"))
# weitere Parameter

# größe der Aggregationsfesnter für die CD Auswertung Distanz-Schwellenwert
# "dist" ist abhängig von wins_size je kleiner win_size desto größer der
# Schwellenwert
win_size = 25
#--- Download der Daten
# gui = TRUE ruft die GUI zur Kontrolle auf
if (get_sen){
  out_paths_3 <- sen2r(
    gui = F,
    param_list = "~/edu/geoinfo/data/harz.json",
    tmpdir = envrmt$path_tmp,
  )
}

#--- Einlesen der Daten aus den Verzeichnissen
# RGB stack der beiden Jahre
pred_stack_2019 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20190619",full.names = TRUE))
pred_stack_2020 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20200623",full.names = TRUE))

# Einlesen der Corine Daten
# Für den Download ist ein Konto notwendig https://land.copernicus.eu/pan-european/corine-land-cover
# Daher die Daten manuell herunterladen und in das Verzeichnis kopieren und entpacken
# Dann auskommentierten Tail ausführen
# corine_eu = raster(file.path(envrmt$path_data_lev0,"u2018_clc2018_v2020_20u1_raster100m/DATA/U2018_CLC2018_V2020_20u1.tif"))
# tmp = projectRaster(pred_stack_2019[[1]],crs = crs(corine_eu))
# corine_crop = raster::crop(corine_eu,tmp)
# corine_utm = projectRaster(corine_crop,crs = crs(pred_stack_2019))
# corine = resample(corine_utm,pred_stack_2019[[1]])
# raster::writeRaster(corine,file.path(envrmt$path_data_lev0,"/corine.tif"),overwrite=TRUE)

# ansonsten den Beipieldatensatz laden
corine = raster::raster(file.path(envrmt$path_data_lev0,"corine.tif"))

# Erstellen einer Wald-Maske
# Agro-forestry areas code=22, Broad-leaved forest code=23,
# Coniferous forest code=24, Mixed forest code=25
mask = reclassify(corine,c(-100,22,0,22,26,1,26,500,0))

# Stack-Loop über die Daten
for (pat in c("RGB432B","EVI","MSAVI2","NDVI","SAVI")){
  pred_stack_2019 = raster::stack(pred_stack_2019,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20190619",full.names = TRUE)))
  pred_stack_2020 = raster::stack(pred_stack_2020,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20200623",full.names = TRUE)))
}
# get rid of NA
pred_stack_2019 = reclassify(pred_stack_2019, cbind(NA, 0))
pred_stack_2020 = reclassify(pred_stack_2020, cbind(NA, 0))

# Zuweisen von leserlichen Namen auf die Datenebenen
names(pred_stack_2019) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
saveRDS(pred_stack_2019,paste0(envrmt$path_data,"pred_stack_2019.rds"))
saveRDS(pred_stack_2020,paste0(envrmt$path_data,"pred_stack_2020.rds"))
pred_stack_2019 = readRDS(paste0(envrmt$path_data,"pred_stack_2019.rds"))
pred_stack_2020 = readRDS(paste0(envrmt$path_data,"pred_stack_2020.rds"))
##---- Kmeans Klassifikation der Daten----
## k-means über RStoolbox
# Verändern Sie den Wert von nclasses und vergleichen sie die Ergebnisse
nclasses = 2
prediction_kmeans_2019 = unsuperClass(pred_stack_2019,
                                      nClasses = nclasses,
                                      norm = TRUE,
                                      algorithm = "MacQueen" )
prediction_kmeans_2020 = unsuperClass(pred_stack_2020,
                                      nClasses = nclasses,
                                      norm = TRUE,
                                      algorithm = "MacQueen")

# Visualisierung mit Mapview: mit dem + Zeichen können beliebige Raster und
# Vektor-Datenebnenen kombiniert werden Bei der Ansicht fällt auf dass
# offensichtlich der Clusteralgorithmus keine zielführenden Klassen der
# Einstellung nClasses=2 erzeugt
mapview(mask*prediction_kmeans_2019$map,
        at = seq(0, nclasses, 1), # Anzahl der Legendenklassen
        legend = TRUE,            # Legende zeigen
        alpha.regions = 1,        # layer undurchsichtig
        maxpixels =  1693870) +  #  volle auflösung
  viewRGB(mask * pred_stack_2019, # funktion viewRGB zeigt dreikanalige Raster Bilder
          r=4,g=5,b=6,            # Kanäle aus dem angegebenen Rasterstack
          maxpixels =  1693870) +
  mapview(mask*prediction_kmeans_2020$map, # RSToolbox Rasterobjekte werden unter $map gespeichert
          at = seq(0, nclasses, 1),
          legend = FALSE,
          alpha.regions = 1,
          maxpixels =  1693870) +  # neue Ebene
  viewRGB(mask * pred_stack_2020,
          r=4,g=5,b=6,
          maxpixels =  1693870)
## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Digitalisierung der Trainingsdaten

#---- Digitalisierung der Trainingsdaten ----
# Muß natürlich nur einmal gemacht werden
if (digitize) {
  source(file.path(root_folder, "src/functions/digitize.R"))
} else {
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
}

Überwachte mit Klassifikation Random Forest

### Maximum Likelihood Classification

# ---- Maximum Likelihood Classification ----
# superClass() Funktion aus dem Paket RSToolbox umwandeln der Tabelle in das
# geforderte SpatialdataPoint Objekt Identifikation der Koordinaten
sp::coordinates(trainDat_2019) = ~X+Y
sp::coordinates(trainDat_2020) = ~X+Y
# Zuweisen der korrekten Projektion (hier aus dem zugehörigen Raster-Stack)
crs(trainDat_2019) = crs(pred_stack_2019)
crs(trainDat_2020) = crs(pred_stack_2020)

# superClass method "mlc" berechnet die Statisik und klassifiziert in einem aufruf
prediction_mlc_2019       <- superClass(pred_stack_2019, trainData = trainDat_2019[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5)
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
prediction_mlc_2020       <- superClass(pred_stack_2020, trainData = trainDat_2020[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5)
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
saveRDS(prediction_mlc_2019,paste0(envrmt$path_data,"pred_mlc_2019.rds"))
saveRDS(prediction_mlc_2020,paste0(envrmt$path_data,"pred_mlc_2020.rds"))
prediction_mcl_2019 = readRDS(paste0(envrmt$path_data,"pred_mlc_2019.rds"))
prediction_mlc_2020 = readRDS(paste0(envrmt$path_data,"pred_mlc_2020.rds"))

# Ergebnisse der rf Klassifikation
prediction_mlc_2019
## superClass results
## ************ Validation **************
## $validation
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut       37    11
##   other           5   670
##                                           
##                Accuracy : 0.9779          
##                  95% CI : (0.9643, 0.9873)
##     No Information Rate : 0.9419          
##     P-Value [Acc > NIR] : 2.471e-06       
##                                           
##                   Kappa : 0.8105          
##                                           
##  Mcnemar's Test P-Value : 0.2113          
##                                           
##             Sensitivity : 0.88095         
##             Specificity : 0.98385         
##          Pos Pred Value : 0.77083         
##          Neg Pred Value : 0.99259         
##              Prevalence : 0.05809         
##          Detection Rate : 0.05118         
##    Detection Prevalence : 0.06639         
##       Balanced Accuracy : 0.93240         
##                                           
##        'Positive' Class : clearcut        
##                                           
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : class 
## values     : 1, 2  (min, max)
## attributes :
##  ID    value
##   1 clearcut
##   2    other
prediction_mlc_2020
## superClass results
## ************ Validation **************
## $validation
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut      114    28
##   other           7   766
##                                           
##                Accuracy : 0.9617          
##                  95% CI : (0.9472, 0.9732)
##     No Information Rate : 0.8678          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8448          
##                                           
##  Mcnemar's Test P-Value : 0.0007232       
##                                           
##             Sensitivity : 0.9421          
##             Specificity : 0.9647          
##          Pos Pred Value : 0.8028          
##          Neg Pred Value : 0.9909          
##              Prevalence : 0.1322          
##          Detection Rate : 0.1246          
##    Detection Prevalence : 0.1552          
##       Balanced Accuracy : 0.9534          
##                                           
##        'Positive' Class : clearcut        
##                                           
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : class 
## values     : 1, 2  (min, max)
## attributes :
##  ID    value
##   1 clearcut
##   2    other
## ---- Visualisierung mit mapview ----
mapview::viewRGB(mask*pred_stack_2020, r = 4, g =5, b = 6,maxpixels =  1693870)+
  mapview(mask*prediction_rf_2019 , alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = TRUE) +
  mapview(mask*prediction_rf_2020, alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2019$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2020$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE)
## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Change Detection Auswertung

Differenzbild random forest

## ---- Change Detection Auswertung ----

# Differenzbild random forest
tmap::qtm(prediction_rf_2020 - prediction_rf_2019,title = "kmeans 2019")
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

# ---- CDA Extraktion der Klassennamen ----
categories = levels(prediction_rf_2019)[[1]]$value
categories
## [1] "clearcut" "other"
# ---- Berechnung der Vierfeld Tabelle mit raster Basisfunktion ----
ct = raster::crosstab(prediction_rf_2019,prediction_rf_2020)
rownames(ct) = categories
colnames(ct) = categories
ct %>%
  kbl(caption = "Crosstab 2019 vs 2020",) %>%
  kable_classic(full_width = F)
(#tab:Detection_Auswertung)Crosstab 2019 vs 2020
clearcut other
clearcut 4716 11640
other 41635 1635879 
# add_header_above(c("layer 1 " = 1, "layer 2" = 2))

Berechnung der Kappa Werte

# ---- kappa ----
# Vergleich der Übereinstimmung unterschiedlicher Klassifikationen (hier MaxLike
# und RF) mit Hilfe diverser Kappa Werte
# https://giswerk.org/doku.php?id=r:r-tutorials:calculatekappa
# 2019
kappa_2019 = kstat(prediction_mlc_2019$map,prediction_rf_2019, perCategory = FALSE)
kappa_2019
##                 K      Kloc    Khisto CA       QA       AA        AD       QD
## overall 0.5285764 0.9298423 0.5684582 50 46.68416 1.752676 0.1322416 1.430924
kappa_2019  %>%
  kbl(caption = "Kappa 2019 MLC vs RF",) %>%
  kable_classic(full_width = F)
Table 1: Kappa 2019 MLC vs RF
K Kloc Khisto CA QA AA AD QD
overall 0.5285764 0.9298423 0.5684582 50 46.68416 1.752676 0.1322416 1.430924 
# 2020
kappa_2020 = kstat(prediction_mlc_2020$map,prediction_rf_2020,perCategory = FALSE)
kappa_2020
##                 K      Kloc    Khisto CA      QA       AA       AD       QD
## overall 0.6513909 0.8387926 0.7765815 50 43.2518 4.395713 0.844811 1.507672
kappa_2020  %>%
  kbl(caption = "Kappa 2019 MLC vs RF",) %>%
  kable_classic(full_width = F) %>%
  column_spec(1,color = "red", link = "https://giswerk.org/doku.php?id=r:r-tutorials:calculatekappa")
Table 1: Kappa 2019 MLC vs RF
K Kloc Khisto CA QA AA AD QD
overall 0.6513909 0.8387926 0.7765815 50 43.2518 4.395713 0.844811 1.507672

Veränderungs Raster aus der Kreuztabelle

# Erzeugen aller Vergleichsraster der Kontingenztabelle
# Die lapply funktionen sind integrierte FOR Schleifen die über die Liste der
# Kategorien die Funktion changefrom() für die Kreuztabelle anwenden
r = lapply(1:length(categories),
           function(i){lapply(1:length(categories),
                              function(j){changefrom(prediction_rf_2019, prediction_rf_2020, i,j)})})
r
## [[1]]
## [[1]][[1]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## [[1]][[2]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## 
## [[2]]
## [[2]][[1]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## [[2]][[2]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
# ---- Visualsierung der Kreuztabellierten Von-Zu-Raster ----
# Plotten der Raster hierzu werden erneut alle Kategorien einzeln geplottet i
# und j sind hilfsvariablen um die korrekten Raster Layer ansprechen zu können.
# t ist eine Hilfvariable um eine Liste für die Ergebnisbilder hochzählen zu
# können
tmap_mode("view")
## tmap mode set to interactive viewing
t=i=j=1 # setze zählvariablen auf 1
m=list() # erzeuge leere Liste für die Karten
for(cat1 in categories)  { # für jede Kategorie
  for(cat2 in categories)  { # mit jeder Kategorie
    # plotte die Karte
    m[[t]]  = tm_shape(st_as_stars(r[[i]][[j]])) +
      tm_raster(col = "layer",palette = "viridis",style = "cat",
                title = if(cat1==cat2) {
                  paste("no changes " , unique(cat1,cat2))
                }
                else if (cat1!=cat2) {
                  paste(cat1," -> ",cat2)
                })
    # zähle die innere schleife hoch
    j = j + 1
    # zähle die ergebnisliste hoch
    t = t + 1
  }
  # innere schleife abgearbeitet setze sie auf  den Anfang zurück
  j = 1
  # zähle die äußere Schleife hoch
  i = i + 1
}
# Interaktive und synchronisierte Karten
tmap::tmap_arrange(m,sync = TRUE)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

Räumliche Analyse der Veränderungen mit dem paket motif

# ---- Analyse von Veränderungen mit dem paket motif ----
# lsp_compare vergleicht zwei (oder mehr) kategoriale Karten (change detection
# klassifikationen etc.) miteinander und nutzt für die Ausgabe der
# Wahrscheinlichkeiten verschiedener räumlicher Aggregierungsstufen und
# Merkmalsraum-Distanzen um Veränderungswahrscheinlichkeiten zu ermitteln
mrf_compare_2020_2019 = lsp_compare(st_as_stars(mask*prediction_rf_2019),
                                    st_as_stars(mask*prediction_rf_2020),
                                    type = "cove",
                                    dist_fun = "jensen-shannon",
                                    window = win_size,
                                    threshold = 0.9)
## Metric: 'jensen-shannon' using unit: 'log2'.
# Visualisierung der Gesamtwahrscheinlichkeiten Achtung logarithmische Skala

tmap_mode("plot")
## tmap mode set to plotting
tm_compare2 = tm_shape(mrf_compare_2020_2019) +
  tm_raster("dist",
            breaks = c(0, 0.001, 0.01, 0.1, 1.01),
            legend.is.portrait = F,
            palette = "viridis",
            title =  "Distance (JSD)") +
  tm_layout(legend.show = TRUE,
            legend.text.size = 0.3,
            legend.outside = TRUE) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = paste("RF 2019 vs 2020 ",win_size*10,"m**2" ),
            title.size = 1.0) +
  tm_grid()

tm_compare2

# ---- Detail-Analyse Teil 1 ----
# Identifikation der Gebiete (Referenz ist win_size) die maximale Veränderungen
# aufweisen im Beispiel soll  "dist" soll größer 0.001 sein (logarithmische
# Skalierung!)
lc_am_compare_sel = st_as_sf(mrf_compare_2020_2019) %>%
  subset(dist > 0.01)

# Sortierung nach Größe
lc_am_compare_sel = lc_am_compare_sel[order(lc_am_compare_sel$dist,decreasing = TRUE), ]
lc_am_compare_sel
## Simple feature collection with 552 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 582400 ymin: 5744920 xmax: 597400 ymax: 5755670
## Projected CRS: WGS 84 / UTM zone 32N
## First 10 features:
##        id na_prop_x na_prop_y      dist                       geometry
## 2107 2107         0         0 0.6196571 POLYGON ((583900 5747420, 5...
## 2291 2291         0         0 0.5594878 POLYGON ((584900 5746670, 5...
## 2290 2290         0         0 0.5212548 POLYGON ((584650 5746670, 5...
## 2470 2470         0         0 0.3918801 POLYGON ((584650 5745920, 5...
## 2471 2471         0         0 0.3341960 POLYGON ((584900 5745920, 5...
## 2650 2650         0         0 0.3203009 POLYGON ((584650 5745170, 5...
## 2231 2231         0         0 0.3052788 POLYGON ((584900 5746920, 5...
## 2531 2531         0         0 0.2897703 POLYGON ((584900 5745670, 5...
## 2116 2116         0         0 0.2495521 POLYGON ((586150 5747420, 5...
## 2649 2649         0         0 0.2460312 POLYGON ((584400 5745170, 5...
##- Visualsierung
# Plotten der top ten gebiete wir nutzen die selbst geschriebene Funktion tm_plot()
# siehe src/functions/fun_tmap.R
tm_plot(sel=lc_am_compare_sel[1:nrow(lc_am_compare_sel),],ov = T)
## tmap mode set to plotting
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.4. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

# als interaktive karte
tm_plot(sel=lc_am_compare_sel[1:nrow(lc_am_compare_sel),],vt = "view",ov = T)
## tmap mode set to interactive viewing
# ---- Visualsierng Detail-Analyse ----
# Extraktion nach der sortierten Liste lc_comapare_sel  werden die top 10 change
# detection hotspots Daten extrahiert und in die liste lc geschrieben
lc = lapply(seq(1:10),function(i){
              lsp_extract(c(st_as_stars(prediction_rf_2019),
                            st_as_stars(prediction_rf_2020)),
                          window = win_size,
                          id = lc_am_compare_sel$id[i])
            })

# Erzeugen der top 3 Change Detection Hotspots Karten-Ausschnitte
tm_lc = lapply(seq(1:3),function(i){
               tm_plot(lc[[i]])
             })
## tmap mode set to plotting
## tmap mode set to plotting
## tmap mode set to plotting
# plotten der erzeugten Karten
tmap_arrange(tm_lc)

#####

Questions and mistakes but also suggestions and solutions are welcome.

Due to an occasionally faulty page redirection, a 404 error may occur. please use the alternative

Content
    Tag Nav
    basics ecology forest forestry lidar project-management remote-sensing reproducibility resources rmarkdown scidoc scripts